Isostaticity of Constraints in Jammed Systems of Soft Prictionless Platonic Solids 
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The average number of constraints per particle (Ctotai) in mechanically stable systems of Platonic 
solids (except cubes) approaches the isostatic limit at the jamming point ({Ctotai) — ► 12), though 
average number of contacts are hypostatic. By introducing angular alignment metrics to classify 
the degree of constraint imposed by each contact, constraints are shown to arise as a direct result 
of local orientational order reflected in edge-face and face-face alignment angle distributions. With 
approximately one face-face contact per particle at jamming chain-like face-face clusters with finite 
extent form in these systems. 

PACS numbers: 45.70.Cc, 45.50.-j, 61.43.-j 



Following Maxwell's approach jammed assemblies 
of frictionlcss spheres exhibit average number of contacts 
per particle (Ztotai) equal to the isostatic value 2n/ 
where n/ is the degrees of freedom per particle. In con- 
trast, ellipses 0,0], ellipsoids 3, tetrahedra @, @], and 
the remaining Platonic solids fj exhibit hypostatic be- 
havior ((Ztotai) < 2nf). The isostatic condition has 
been linked to the mechanical stability of soft sphere 
systems 0], and the hypostaticity of ellipses has been 
attributed to the presence of floppy vibrational modes, 
which provide vanishing restoring force Q. Jaoshvili 
et al. [5[ recently asserted that the average constraint 
number (Ctotai), which incorporates topology-dependent 
contact constraint (e.g., through vertex-face, edge-edge, 
edge- face, and face-face contact topologies), is isostatic 
for tetrahedral dice even though (Z to tai) < 2n/. Simi- 
lar approaches incorporating variable contact constraint 
have been utilized in the prediction of mechanism mobil- 
ity as early as 1890 (see [8j). The presupposition under- 
lying constraint counting as a valid means of predicting 
net degrees of freedom is the mechanical independence 
of constraints. As a result, such constraint counting ap- 
proaches fail to predict the true mobility of mechanisms 
when kinematic redundancy is present (see 0), and iso- 
static (Ctotai) is not guarranteed at the jamming point. 
Though constraint counting may yield an isostatic result, 
the validity of this conclusion depends on the methods 
used to estimate (Ctotai) as we show in this Letter. 

The main goal of this Letter is to assess the validity 
of generalized isostaticity Q for jammed systems of Pla- 
tonic solids (except cubes) through energy-based struc- 
tural optimization and objective topological classifica- 
tion. We also present jamming threshold density and 
show that these packings are vibrationally stable. The 
ill-conditioned vibrational spectrum is a direct result of 
topological heterogeneity in the contact network, and we 



therefore identify and highlight the orientational order 
of contacts by introducing angular alignment metrics. 
Topology classification via angular alignment metric dis- 
tributions is thereby used to assess the isostaticity of con- 
straints. 

Contact model and jamming protocol. Structural opti- 
mization coupled with controlled consolidative and ex- 
pansive strain is used to probe the jamming point as 
in The conservative model employed assumes that 
contact between particles a and (3 results in energy 
E a p = YV 2 /AV P after a Hookian contact model applied 
to uniaxially compressed bars, where V is the intersec- 
tion volume between the particles, V p is the volume of a 
single particle, and Y is the elastic modulus. Conjugate 
gradient minimization with line searching is utilized with 
a relative energy change convergence tolerance less than 
10 -12 at each strain step to simulate static equilibrium 
(see [(§). Density is defined as <f> — NV p /V ce u, where V ce ii 
and N are the volume and the number of particles in the 
primary periodic cell. 



TABLE I: Jamming threshold 
with 95 % confidence intervals. 



)j estimated for each system 



A' = 100 



N = 400 



exp. [9J° 



Tetrahedra 0.629±0.001 0.634±0.011 0.64±0.01 

Octahedra 0.6796±0.0003 0.686±0.001 0.64±0.01 

Icosahedra 0.6953±0.0003 0.7008±0.0003 0.59±0.01 

Dodccahcdra 0.7065±0.0002 0.7085±0.0003 0.63±0.01 

"Values presented in [|| for tetrahedral dice. 

Assemblies of monodisperse Platonic solids with peri- 
odic cubic lattice boundary conditions were consolidated 
with average energy per particle near 3.2 x 10~ 5 YV P from 
low-density random configurations at <j> = 0.05. Esti- 
mates of the jamming threshold density <pj (Table I) were 
obtained by expansion toward the jamming point as in 
4>j converges well at N = 100, but the results differ 
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FIG. 1: (Color online) (a) Vibrational spectra for stable sys- 
tems nearest the jamming threshold, (b) The power law scal- 
ing of (uji) (lower curve set) and the asymptotically constant 
scaling of (uj h ) (upper curve set). 



somewhat from the 'random close packed' densities mea- 
sured by [§| for finite systems of rounded dice. Unless 
otherwise stated, hereafter the results presented are for 
N = 400. 

Mechanical stability. We computed the low-energy 
vibrational spectra of these monodisperse systems as 
LUi = yj Aj Jrn (as in Q), where Xi is the i th eigen- 
value of the dynamical matrix D a p = d 2 E / dr a drp 
and m is particle mass. Frequency uii is presented in 



units of V p y/Y/p, where p is mass density of the 
solid phase. Coordinates r a of particle a are com- 
posed of translational x% t0t and rotational 6^ a com- 
ponents r a = {xi ta ,X2, a ,X3 ta ,R a 6i !a ,R a 62, a ,Ra^3,a}: 

where R a is the radius of gyration. We calculate D a p 
through central differences of forces and moments, con- 
sidering only its symmetric part. The resulting spectra 
of static equilibrium systems with density nearest <f>j are 
displayed in Fig. 1(a). For static equilibrium systems at 
all densities we find 6N — 3 stable modes with 3 triv- 
ial (rigid body) translational modes, confirming that our 
soft packings are indeed stably jammed as in 

UM- This 

result also confirms that all N particles participate in the 
mechanical network, and therefore no rattlers exist. 

The mean values of the 50 lowest and highest non- 
trivial frequency modes (w;) and (ujh) were computed and 
plotted against excess density, Acp = <p — 4>J [Fig- 1(b)]. 
The scaling of (u>i) with respect to A</>, (wj) ~ (A</>) 2 , 
reveals that these packings are in fact marginally stable 
at <j>j, while finite (w/j) persists at the jamming point 
[Fig. 1(b)]. The latter behavior is a signature of trans- 
lational vibrational modes involving face-face contacts, 
because only such contact topologies exhibit harmonic 
(i.e., quadratic) energy variation and consequently con- 
stant stiffness with respect to displacements along the 
direction of contact force induced by straining from <f>j. 
By classifying contact topologies with methods outlined 
subsequently we have confirmed that these modes are lo- 
calized on clusters formed by face-face contacts. Energies 
of these systems scale as E oc (A</>) 6 Q with bulk mod- 
ulus K scaling as K oc (A</>) 4 . Therefore it is clear that 
low frequency modes are excited by volumetric strain be- 
cause (u){) oc \[K oc (A0) 2 . 



Angular alignment. Soft contacts between faceted par- 
ticles can be classified as face-face, edge-face, vertex-face, 
or edge-edge. Hereafter vertex-face and edge-edge con- 
tacts are referred to as 'lower order' since they exhibit less 
order relative to face-face and edge-face contacts. Ac- 
cordingly we define angular order metrics that approach 
zero as contacts orient with perfect face-face or edge- 
face alignment; these metrics could be measured experi- 
mentally with tomographic reconstruction. The face-face 
alignment angle 0f-f.qi of contacting particles q and I 
is expressed as 8f-f >q i — n — cos _1 (min(nJny)), where 
h q i is the normal vector of face i on particle q. The 
minimum is computed over all combinations of i and 
j corresponding to intersecting faces on the respective 
particles. The edge-face alignment angle # e _/ is cal- 
culated as the minimum of 9 e -f^ q i and 6 e —fjq, where 
@e-f,qi = sin - (mia(|e^riy|)), and e q i is the unit edge 
vector of edge i on particle q. 

Randomly oriented faces and edges provide a start- 
ing point for understanding alignment angle distribu- 
tions in jammed systems. Three-dimensional (3D) ran- 
dom edges and faces possess probability density of edge- 
face alignment with p(9 e -f) oc cos(# e _/). Edge-face con- 
tacts arc therefore expected to be ubiquitous in jammed 
systems, because probability is weighted toward small 
6 e -f. In contrast, face-face contacts are expected to be 
less common, because 3D random faces possess align- 
ment probability density that vanishes in the face-face 
limit \p(6f-f) oc sin(0/_/)]. On the other hand, p(6f-f) 
for edge-face constrained contacts exhibit uniform non- 
vanishing probability density. 



TABLE II: Critical edge-face c -[,, ., and face-face 9j -/., ., 
alignment angles for each topology i in degrees. 
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Cumulative distribution functions (CDFs) are plot- 
ted in Fig. 2 for all the contacts in one realization of 
each system in addition to CDFs of ideal random sys- 
tems. Simulated CDFs exhibit critical angles plotted in 
Fig. 2 and listed in Table II that bound the possible 
alignment angles for face-vertex, edge-edge, and edge- 
face contact topologies; we denote critical edge-face and 
face-face alignment angles as # e _/ jCj! ; and 9f-f tC ,i, respec- 
tively, where i denotes the particular contact topology. 
These angles correspond to pairs of particles contacting 
with a particular topology oriented with highest possi- 
ble symmetry. The angular breadth of c(O e -f) spans 
proportional to 9 e ~f tCt mi n , the minimum critical edge- 
face alignment angle among vertex-face and edge-edge 
contacts. This results in particles with large e -/,c,)ni« 
having low probability density (i.e., CDF slope) of near 
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FIG. 2: (Color online) CDFs of (a) 6 e . f and (b) 9 J - J for 
stable systems nearest <j>j. Points plotted on the simulated 
curves correspond to c(0 e _/ lC ,i) and c(f9/_/ iCl j). Region 1 of 
c(Qj-f) is masked in gray, while region 2 extends from the 
edge of region 1 to #/-/,c,roinj the minimum critical face- 
face alignment angle. Ideal random CDFs for edge-face con- 
strained (open squares) and unconstrained (open circles) con- 
tacts are displayed as well. 



cdgc-facc contacts. All systems exhibit hyper-random 
probability density at small 9 e —f with CDF approaching 
that of edge-face constrained contacts [open squares, Fig. 
2(a)]. 

Fig. 2(b) shows that c(0/_/) for all systems possess 
three regions of descending probability density (i.e., CDF 
slope): (1) ultra- low angle (< 1°) face- face contact re- 
gion, (2) intermediate uniform probability region, and 
(3) high angle region. Region 2 contains the largest por- 
tion of total contacts and exhibits c(0/_/) similar to that 
of random edge-face constrained contacts [open squares, 
Fig. 2(b)], rather than unconstrained contacts [open cir- 
cles, Fig. 2(b)]. This results from the abundance of con- 
tacts with small Q e -t that exhibit CDFs approaching 
that of edge-face constrained contacts [open squares, Fig. 
2(a)], rather than unconstrained contacts [open circles, 
Fig. 2(a)]. Region 2 possesses approximately invariant 
p(9f-f) resulting in increased probability for contacts to 
form with acute 0f-f [Fig. 2(b)] relative to completely 
unconstrained random systems [open circles, Fig. 2(b)]. 
The emergence of region 1 [gray mask, Fig. 2(b)] can be 
understood by considering the constrained rotation of an 
ideal edge-face contact. During consolidation such a con- 
tact will rotate randomly about its edge through a range 
of #/-/• Repulsion between opposing faces at #/-/ = 
will prevent the contact from rotating further. We ob- 
serve the evolution of c(0 e —f) and c(#/_/) as <p ~ ► 4>J 
qualitatively consistent with this idealized picture - re- 
gion 1 of c(6f-f) emerges only after c(6 e -f) becomes 
hyper-random. 

Constraint isostaticity and topological distributions. 
To classify topologies we fit c(0 e _/) and c(9f-f) indepen- 
dently with a piecewise continuous distribution function 
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basis, Cf(8) = H(9 — a)[aQ — o n (9 — a) n ], where a n are 

71=1 

fitted coefficients, H is the Heaviside step function, and 
a is fitted CDF discontinuity. All empirical CDFs have 
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FIG. 3: (Color online) Variation of average (a) contact, (b) 
constraint, and (c) face- face contact number with A<j>. Half- 
filled symbols represent data for systems of TV = 100, while 
the remainder are for N = 400. Data points of tetrahedra 
and octahedra at Acj> = 0.01 and tetrahedra at A<f> = 0.02 ex- 
hibited mild vibrational instability (5.3 % unstable modes for 
tetrahedra at A<f> — 0.01 and less than 1.2 % for the others), 
(d) Face-face clusters with maximal linear extent formed in 
stable TV = 400 systems nearest <f>j of octahedra, tetrahedra, 
ieosahedra, and dodecahedra (from left to right). 



where i corresponds to the particular alignment angle. 
With the fitted parameter do angular alignment cutoffs 
for topological classification 9 cut are determined such 
that c(9 cu t) = ao- Edge- face contacts are classified as 



those with 9 e -f < 9 e -f,cut an d 0/-/ > 
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face contacts as those with 0/-/ < 0/-/,c«t- 

The variation of average contact number (Z total) with 
respect to A</> [Fig. 3(a)] confirms the generally hypo- 
static nature of (Z to tai), consistent with our previous 
findings for smaller systems @. We determine the av- 
erage contact number for face-face (Zf-f), edge- face 
(Z e —f), and lower order contacts (Zi). The average con- 
straint number (C total) is thereby calculated as (C 'total) = 
3{Zf- f ) + 2(Z e -f) + (Zi) 0. Importantly, (C to tai) of 
each system approaches values near the isostatic limit 
[Fig. 3(b)], in contrast with (Z to tai) [Fig. 3(a)]. Thus, 
these systems possess contacts that independently con- 
strain motion. These systems also possess (Zf-f) < 1 
[Fig. 3(c)]. We attribute this effect to a two-fold rota- 
tional constraint induced on a given particle once a single 
face-face contact is formed. Such a rotational constraint 
appears to hinder the formation of additional face-face 
contacts. 

Note that the low (Zf-f) of tetrahedra clearly con- 
trasts with the value of 2.3 recently reported for tetrahe- 
dral dice The practical importance of this finding is 
very significant, for if two to three face- face contacts per 
particle were present, as reported by we expect such 
systems to readily exhibit cluster percolation and radi- 
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TABLE III: Distributions of average contact and constraint 
numbers for the present work and prior work. 
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11.1 


6.3 


(Zf-f) 


0.9 
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2.3 


(Ze-f) 


1.5 


3.7 


1.2 


(z,) 


8.6 


2.0 


2.8 


(C total ) 


14.5 


25.5 


12.1 



a 100 tetrahedra at A0 = 0.04. 

''Values are computed by averaging over results obtained with 
6Ut = {5,15,25,35,45,55} as in gEfl. 
c Values presented in [5|] for tetrahedral dice. 

cally different structures and mechanical behavior than 
those with near unity {Zf-f). In Table III we compare 
the contact and constraint numbers of the present work 
(column 1) with that of Q (column 3). Indeed, all con- 
tact numbers are markedly different (except (Z e -f)) from 
that of ■ We also utilize a topological classification 
procedure similar to that in [H], [jj by averaging contact 
numbers obtained for various values of 6 cut (column 2 
in Table III) . The resulting constraint number is twice as 
large as the isostatic value. Thus, arbitrary choice of cut , 
as employed in [B], [l^ , can yield a geometrically infeasi- 
ble range of contact numbers and generally non-isostatic 

(Ctotal)- 

TABLE IV: Face-face contact number (Zf-f). average cluster 
size S, maximal extent along Cartesian axes Imax, and the 
fractal dimension of clusters with maximal linear extent D„ lax 
for stable N = 400 systems nearest <f)j. 

Octahcdra Icosahcdra Tetrahedra Dodccahcdra 

(Zf-f) 059 077 5794 1702 

S 1.88 2.45 3.18 4.13 

i m a*/(V„n) 1/3 0.665 0.795 0.741 0.776 

D mal 1.01±0.10 1.23±0.08 1.37±0.09 1.32±0.14 



from percolation. Thus, the lack of percolation is a di- 
rect consequence of near unity (Zf-f). These results 
suggest that face-face cluster formation is a bond per- 
colation process with respect to (Zf-f). Therefore, we 
expect that systems of particles with shapes conducive 
to ordering (e.g., cubes) or with attractive interactions 
could increase (Zf-f), forming larger clusters that per- 
colate. 

In summary, we have established the isostaticity of 
constraints in disordered jammed systems of all Platonic 
solids except cubes and have linked this condition to their 
mechanical stability. Our results suggest that (Zf-f) or 
other integral functions of p(9f-f) are suitable order pa- 
rameters that can be used to determine the maximally 
random jammed state of faceted particle systems accord- 
ing to the approach described in []]| . The structure and 
extent of face-face clusters is found to be a consequence 
of few face- face contacts in these systems. Future work 
will focus on identifying means by which face-face contact 
number may increase and studying the critical behavior 
of such systems with ensuing face-face contact percola- 
tion. 
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We have analyzed the structure of clusters formed by 
face- face contacts (Table IV). The topological connec- 
tivity of clusters was considered under periodic bound- 
ary conditions. The average size of clusters S is defined 
in terms of the number of particles s in each cluster as 
S = J2 s 2 / s > we nn d that S increases with (Zf-f) as 
particle shape is changed. We also find that all clusters 
in systems of N = 400 particles do not percolate ac- 
cording to topological connectivity. Such a requirement 
for periodic percolating networks is stricter than the re- 
quirement that a cluster spans the system boundaries 
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Therefore, we have also considered the less restric- 
tive percolation requirement that a cluster spans system 
boundaries. Clusters in the respective systems exhibit 
maximal extent along Cartesian axes l m ax smaller than 
the length of the finite cubic cell. Clusters with maximal 
linear extent exhibit chain-like structures [Fig. 3(d)] with 
fractal dimension D max > 1 (Table IV). Considering the 
fractal dimension of ~ 2.5 for a percolating cluster at 
the threshold in a simple cubic lattice [13], these chains 
are substantially lower dimensional and appear to be far 
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